
*****************************************************************************
****************************** 24 Month Lags ********************************
*****************************************************************************

clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO


** Controls in LP regression **

* 24 months of lags of Federal Funds 
local MaxLPLags = 24
* 12 months of lags of dependent variable
local MaxINVLags = 24


* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').FF l(1/`MaxINVLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').FF l(1/`MaxINVLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').FF l(1/`MaxINVLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').FF l(1/`MaxINVLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').FF l(1/`MaxINVLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.03 38 "Abstract", color(green) size(small)) ///
	text(-.03 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI24, replace) legend(off)
	graph save "${path}\Results\Robust\WLI24.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.002 38 "Abstract", color(green) size(small)) ///
	text(-.007 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(RW24, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage24.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.026 38 "Abstract", color(green) size(small)) ///
	text(-.015 38 "Routine", color(navy) size(small)) text(-.004 38 "Manual", color(black) size(small)) ///
	name(EMP24, replace) legend(off)
	graph save "${path}\Results\Robust\Employment24.gph", replace
	
*


******************************************************************************
***************************** 36 Months Lags *********************************
******************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 36

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').FF l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').FF l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').FF l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').FF l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').FF l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI36, replace) legend(off)
	graph save "${path}\Results\Robust\WLI36.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RW36, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage36.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMP36, replace) legend(off)
	graph save "${path}\Results\Robust\Employment36.gph", replace
	
*



*****************************************************************************
******************************** Post 1984 **********************************
*****************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO

keep if year >= 1984


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 12

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').FF l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').FF l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').FF l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').FF l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').FF l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLIpost84, replace) legend(off)
	graph save "${path}\Results\Robust\WLI_post84.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RWpost84, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage_post84.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMPpost84, replace) legend(off)
	graph save "${path}\Results\Robust\Employment_post84.gph", replace
	
*


*****************************************************************************
*************************** Post 1984 w/ 24 lags ****************************
*****************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO

keep if year >= 1984


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 24

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').FF l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').FF l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').FF l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').FF l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').FF l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI24post84, replace) legend(off)
	graph save "${path}\Results\Robust\WLI24_post84.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RW24_post84, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage24_post84.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMP24post84, replace) legend(off)
	graph save "${path}\Results\Robust\Employment24_post84.gph", replace
	
*

*****************************************************************************
*************************** Post 1984 w/ 36 lags ****************************
*****************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO

keep if year >= 1984


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 36

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').FF l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').FF l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').FF l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').FF l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').FF l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI36post84, replace) legend(off)
	graph save "${path}\Results\Robust\WLI36_post84.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RW36post84, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage36_post84.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMP36post84, replace) legend(off)
	graph save "${path}\Results\Robust\Employment36_post84.gph", replace
	
*



*****************************************************************************
*************************** Shocks as Controls ******************************
*****************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO


** Controls in LP regression **

* 24 months of lags of Federal Funds 
local MaxLPLags = 12
* 12 months of lags of dependent variable
local MaxINVLags = 12


* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').RR l(1/`MaxINVLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').RR l(1/`MaxINVLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').RR l(1/`MaxINVLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').RR l(1/`MaxINVLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').RR l(1/`MaxINVLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.03 38 "Abstract", color(green) size(small)) ///
	text(-.03 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\WLI12_shocks.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.002 38 "Abstract", color(green) size(small)) ///
	text(-.007 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(RW_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage12_shocks.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.026 38 "Abstract", color(green) size(small)) ///
	text(-.015 38 "Routine", color(navy) size(small)) text(-.004 38 "Manual", color(black) size(small)) ///
	name(EMP_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\Employment12_shocks.gph", replace
	
*


*****************************************************************************
****************************** 24 Month Lags ********************************
*****************************************************************************

clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO


** Controls in LP regression **

* 24 months of lags of Federal Funds 
local MaxLPLags = 24
* 12 months of lags of dependent variable
local MaxINVLags = 24


* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').RR l(1/`MaxINVLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').RR l(1/`MaxINVLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').RR l(1/`MaxINVLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').RR l(1/`MaxINVLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').RR l(1/`MaxINVLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.03 38 "Abstract", color(green) size(small)) ///
	text(-.03 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI24_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\WLI24_shocks.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.002 38 "Abstract", color(green) size(small)) ///
	text(-.007 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(RW24_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage24_shocks.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.026 38 "Abstract", color(green) size(small)) ///
	text(-.015 38 "Routine", color(navy) size(small)) text(-.004 38 "Manual", color(black) size(small)) ///
	name(EMP24_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\Employment24_shocks.gph", replace
	
*


******************************************************************************
***************************** 36 Months Lags *********************************
******************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 36

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').RR l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').RR l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').RR l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').RR l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').RR l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI36_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\WLI36_shocks.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RW36_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage36_shocks.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMP36_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\Employment36_shocks.gph", replace
	
*



*****************************************************************************
******************************** Post 1984 **********************************
*****************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO

keep if year >= 1984


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 12

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').RR l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').RR l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').RR l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').RR l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').RR l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLIpost84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\WLI_post84_shocks.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RWpost84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage_post84_shocks.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMPpost84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\Employment_post84_shocks.gph", replace
	
*


*****************************************************************************
*************************** Post 1984 w/ 24 lags ****************************
*****************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO

keep if year >= 1984


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 24

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').RR l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').RR l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').RR l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').RR l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').RR l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI24post84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\WLI24_post84_shocks.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RW24post84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage24_post84_shocks.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMP24post84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\Employment24_post84_shocks.gph", replace
	
*

*****************************************************************************
*************************** Post 1984 w/ 36 lags ****************************
*****************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO

keep if year >= 1984


** Controls in LP regression **

* 12 months of lags
local MaxLPLags = 36

* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen RW`i'    = f`i'.RW - l1.RW
	gen WLI`i'   = f`i'.WLI - l1.WLI
	gen Emp`i'   = f`i'.Emp - l1.Emp
	gen hours`i' = f`i'.hours-l1.hours
	gen Output`i' = f`i'.Output-l1.Output
	
}

forvalues v = 1/`o' {
gen b_`v'WLI=.
gen b_`v'RW=.
gen b_`v'Emp=.
gen b_`v'hours=.


gen se_`v'WLI=.
gen se_`v'RW=.
gen se_`v'Emp=.
gen se_`v'hours=.

}

gen b_Output=.
gen se_Output=.

local rhsWeeklyLaborIncome l(1/`MaxLPLags').RR l(1/`MaxLPLags').WLI time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
						
local rhsRealWage l(1/`MaxLPLags').RR l(1/`MaxLPLags').RW time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
					
local rhsEmp l(1/`MaxLPLags').RR l(1/`MaxLPLags').Emp time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 
	
local rhshours l(1/`MaxLPLags').RR l(1/`MaxLPLags').hours time time2
*l(1/`MaxLPLags').logIP l(1/`MaxLPLags').logCPI 

local rhsOutput l(1/`MaxLPLags').RR l(1/`MaxLPLags').Output time time2

**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RR `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLI  = _b[RR] if _n == `i'+1
	replace se_`v'WLI = _se[RR] if _n == `i'+1
			
	* Real Wage LP
	newey RW`i' RR  `rhsRealWage' if OCC == `v', lag(3)
	
	replace b_`v'RW  = _b[RR] if _n == `i'+1
	replace se_`v'RW = _se[RR] if _n == `i'+1
	
	* Employment LP 
	newey Emp`i' RR `rhsEmp' if OCC == `v', lag(3)
				
	replace b_`v'Emp  = _b[RR] if _n == `i'+1
	replace se_`v'Emp = _se[RR] if _n == `i'+1
	
	* Hours LP 
	newey hours`i' RR `rhshours' if OCC == `v', lag(3)
				
	replace b_`v'hours  = _b[RR] if _n == `i'+1
	replace se_`v'hours = _se[RR] if _n == `i'+1
	
			}
}

forvalues i=0/`horizon' {	
	
	* LP regression

	* Output LP 
	newey Output`i' RR `rhsOutput' if OCC == 1, lag(3)
				
	replace b_Output  = _b[RR] if _n == `i'+1
	replace se_Output = _se[RR] if _n == `i'+1
	
			}


****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLI = b_`v'WLI + invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn_`v'WLI = b_`v'WLI - invnormal(1-sig1/2)*se_`v'WLI if _n <= (`horizon' + 1)

	gen up2_`v'WLI = b_`v'WLI + invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	gen dn2_`v'WLI = b_`v'WLI - invnormal(1-sig2/2)*se_`v'WLI if _n <= (`horizon' + 1)
	
	gen up_`v'RW = b_`v'RW + invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn_`v'RW = b_`v'RW - invnormal(1-sig1/2)*se_`v'RW if _n <= (`horizon' + 1)

	gen up2_`v'RW = b_`v'RW + invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	gen dn2_`v'RW = b_`v'RW - invnormal(1-sig2/2)*se_`v'RW if _n <= (`horizon' + 1)
	
	gen up_`v'Emp = b_`v'Emp + invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn_`v'Emp = b_`v'Emp - invnormal(1-sig1/2)*se_`v'Emp if _n <= (`horizon' + 1)

	gen up2_`v'Emp = b_`v'Emp + invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	gen dn2_`v'Emp = b_`v'Emp - invnormal(1-sig2/2)*se_`v'Emp if _n <= (`horizon' + 1)
	
	gen up_`v'hours = b_`v'hours + invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn_`v'hours = b_`v'hours - invnormal(1-sig1/2)*se_`v'hours if _n <= (`horizon' + 1)

	gen up2_`v'hours = b_`v'hours + invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)
	gen dn2_`v'hours = b_`v'hours - invnormal(1-sig2/2)*se_`v'hours if _n <= (`horizon' + 1)

	}

	gen up2_Output = b_Output + invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)
	gen dn2_Output = b_Output - invnormal(1-sig2/2)*se_Output if _n <= (`horizon' + 1)

	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI RW Emp hours {
	forvalues v = 1/`o' {
	
	replace b_`v'`var' = -.25*b_`v'`var'
	replace up_`v'`var' = -.25*up_`v'`var'
	replace dn_`v'`var' = -.25*dn_`v'`var'
	replace up2_`v'`var' = -.25*up2_`v'`var'
	replace dn2_`v'`var' = -.25*dn2_`v'`var'
	
	}
}

*

save "${path}/Clean Data/Impulse Responses.dta", replace

/************************ Impulse Response Functions *********************************/


**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLI Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLI Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLI Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLI Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLI Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLI Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.015 38 "Abstract", color(green) size(small)) ///
	text(-.012 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLI36_post84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\WLI36_post84_shocks.gph", replace
	
************
* Real Wage	
************

*
twoway 	(line b_1RW Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2RW Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3RW Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3RW Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1RW Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2RW Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Average Real Wage", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.004 38 "Abstract", color(green) size(small)) ///
	text(-.005 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(RW36post84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\RealWage36_post84_shocks.gph", replace
	

  ***********
  *Employment
  ***********
  
  *

twoway 	(line b_1Emp Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2Emp Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3Emp Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) ///
	(line up2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3Emp Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1Emp Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2Emp Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.03(.015).03) ytitle("Log Total Employment", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.014 38 "Abstract", color(green) size(small)) ///
	text(-.008 38 "Routine", color(navy) size(small)) text(.001 38 "Manual", color(black) size(small)) ///
	name(EMP36post84_shocks, replace) legend(off)
	graph save "${path}\Results\Robust\Employment36_post84_shocks.gph", replace
	
*

***********************************************************************************
***************************** Asymmetric Effects **********************************
***********************************************************************************


clear 
set more off

use "${path}\Clean Data\CPS Final.dta"

sum FF, detail

* Mean FF is 6.5 percent * 

gen RRhigh = 0
replace RRhigh = RR if FF > 6.5

gen RRlow = 0
replace RRlow = RR if FF <= 6.5


rename WeeklyLaborIncome WLI
rename RealWage RW
rename Emp EMP

xtset OCC time

gen time2 = time*time

* Decide to use original shocks, or shocks estimated with GARCH
drop RR
rename resid_garch RR

egen t = min(time)

replace RW = log(RW)
replace WLI = log(WLI)
gen Emp = log(emp)
gen hours = log(Hours)

gen Output = INDPRO


** Controls in LP regression **

* 24 months of lags of Federal Funds 
local MaxLPLags = 12
* 12 months of lags of dependent variable
local MaxINVLags = 12


* 3 year horizon
local horizon 36

* Romer and Romer shocks 
local policy RR

* O is the number of occupation subgroups under consideration
local o 3

***** LHS variable: the response variable

forvalues i=0/`horizon' {
	
	gen WLI`i'   = f`i'.WLI - l1.WLI
	
}

forvalues v = 1/`o' {
	
gen b_`v'WLIhigh=.
gen b_`v'WLIlow=.

gen se_`v'WLIhigh=.
gen se_`v'WLIlow=.

}

local rhsWeeklyLaborIncome l(1/`MaxLPLags').RR l(1/`MaxINVLags').WLI time time2
						
**
					    
****************************************************************
*** LP table Using Basic OLS
****************************************************************
						
forvalues v = 1/`o' { 
forvalues i=0/`horizon' {	
	
	* LP regression

	* WLI LP 
	newey WLI`i' RRhigh RRlow `rhsWeeklyLaborIncome' if OCC == `v', lag(3)
				
	replace b_`v'WLIhigh  = _b[RRhigh] if _n == `i'+1
	replace se_`v'WLIhigh = _se[RRhigh] if _n == `i'+1
	
	replace b_`v'WLIlow  = _b[RRlow] if _n == `i'+1
	replace se_`v'WLIlow = _se[RRlow] if _n == `i'+1
			
			}
}



****************************************************************
*** Baseline OLS LP graphs
****************************************************************

gen Months = _n-1 if _n <= `horizon' +1

* zero line
gen zero = 0 if _n <= `horizon' +1


***** create confidence bands (in this case 90 and 95%) ****
	scalar sig1 = 0.05	 // specify significance level
	scalar sig2 = 0.3	 // specify significance level
	
	forvalues v = 1/`o' {

	gen up_`v'WLIhigh = b_`v'WLIhigh + invnormal(1-sig1/2)*se_`v'WLIhigh if _n <= (`horizon' + 1)
	gen dn_`v'WLIhigh = b_`v'WLIhigh - invnormal(1-sig1/2)*se_`v'WLIhigh if _n <= (`horizon' + 1)

	gen up2_`v'WLIhigh = b_`v'WLIhigh + invnormal(1-sig2/2)*se_`v'WLIhigh if _n <= (`horizon' + 1)
	gen dn2_`v'WLIhigh = b_`v'WLIhigh - invnormal(1-sig2/2)*se_`v'WLIhigh if _n <= (`horizon' + 1)
	
	gen up_`v'WLIlow = b_`v'WLIlow + invnormal(1-sig1/2)*se_`v'WLIlow if _n <= (`horizon' + 1)
	gen dn_`v'WLIlow = b_`v'WLIlow - invnormal(1-sig1/2)*se_`v'WLIlow if _n <= (`horizon' + 1)

	gen up2_`v'WLIlow = b_`v'WLIlow + invnormal(1-sig2/2)*se_`v'WLIlow if _n <= (`horizon' + 1)
	gen dn2_`v'WLIlow = b_`v'WLIlow - invnormal(1-sig2/2)*se_`v'WLIlow if _n <= (`horizon' + 1)

	}
	

* Flip the impulse response functions to have a negative shock 

foreach var in WLI {
	forvalues v = 1/`o' {
	
	replace b_`v'`var'high = -.25*b_`v'`var'high
	replace up_`v'`var'high = -.25*up_`v'`var'high
	replace dn_`v'`var'high = -.25*dn_`v'`var'high
	replace up2_`v'`var'high = -.25*up2_`v'`var'high
	replace dn2_`v'`var'high = -.25*dn2_`v'`var'high
	
	replace b_`v'`var'low = -.25*b_`v'`var'low
	replace up_`v'`var'low = -.25*up_`v'`var'low
	replace dn_`v'`var'low = -.25*dn_`v'`var'low
	replace up2_`v'`var'low = -.25*up2_`v'`var'low
	replace dn2_`v'`var'low = -.25*dn2_`v'`var'low
	
	}
}

*

**********************
* Weekly Labor Income
**********************

twoway 	(line b_1WLIhigh Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLIhigh Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLIhigh Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLIhigh Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLIhigh Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLIhigh Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLIhigh Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLIhigh Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLIhigh Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.045(.015).045) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.001 38 "Abstract", color(green) size(small)) ///
	text(.011 38 "Routine", color(navy) size(small)) text(.006 38 "Manual", color(black) size(small)) ///
	name(WLIhigh, replace) legend(off)
	graph save "${path}\Results\Robust\WLIhigh.gph", replace


twoway 	(line b_1WLIlow Months, lcolor(black) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_2WLIlow Months, lcolor(navy) ///
	lpattern(solid) lwidth(medium)) /// 
	(line b_3WLIlow Months, lcolor(green) ///
	lpattern(solid) lwidth(medium)) /// 
	(line up2_3WLIlow Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_3WLIlow Months, lcolor(eltgreen) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_1WLIlow Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_1WLIlow Months, lcolor(gray) ///
	lpattern(dash) lwidth(medium)) ///
	(line up2_2WLIlow Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)) ///
	(line dn2_2WLIlow Months, lcolor(eltblue) ///
	lpattern(dash) lwidth(medium)), ///
	ylabel(-0.045(.015).045) ytitle("Log Total Weekly Labor Income", size(medsmall)) xtitle("Months", size(medsmall)) ///
	graphregion(color(white)) plotregion(color(white)) text(.016 38 "Abstract", color(green) size(small)) ///
	text(-.017 38 "Routine", color(navy) size(small)) text(.004 38 "Manual", color(black) size(small)) ///
	name(WLIlow, replace) legend(off)
	graph save "${path}\Results\Robust\WLIlow.gph", replace



	